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We combine the Density Matrix Technique (DMRG) with Green Function Monte Carlo (GFMC) 
simulations. Both methods aim to determine the groundstate of a quantum system but have diflterent 
limitations. The DMRG is most successful in 1-dimensional systems and can only be extended to 
2-dimensional systems for strips of limited width. GFMC is not restricted to low dimensions but 
is limited by the efficiency of the sampling. This limitation is crucial when the system exhibits a 
so-called sign problem, which on the other hand is not a particular obstacle for the DMRG. We 
show how to combine the virtues of both methods by using a DMRG wavefunction as guiding wave 
function for the GFMC. This requires a special representation of the DMRG wavefunction to make 
the simulations possible within reasonable computational time. As a test case we apply the method 
to the 2-dimensional frustrated Heisenberg antiferromagnet. By supplementing the branching in 
GFMC with Stochastic Reconfiguration (SR) we get a stable simulation with a small variance also 
in the region where the ffuctuations due to minus sign problem are maximal. The sensitivity of the 
results to the choice of the guiding wavefunction is extensively investigated. 

We analyse the model as a function of the ratio of the next -nearest to nearest neighbor coupling 
strength which is a measure for the frustration. In agreement with earlier calculations it is found 
from the DMRG wavefunction that for small ratios the system orders as a Neel type antiferromagnet 
and for large ratios as a columnar antiferromagnet. The spin stiffness suggests an intermediate 
regime without magnetic long range order. The energy curve indicates that the columnar phase is 
separated from the intermediate phase by a first order transition. The combination of DMRG and 
GFMC allows to substantiate this picture by calculating also the spin correlations in the system. We 
observe a pattern of the spin correlations in the intermediate regime which is in-between dimerlike 
and plaquette type ordering, states that have recently been suggested. It is a state with strong 
dimerization in one direction and weaker dimerization in the perpendicular direction and thus it 
lacks the the square symmetry of the plaquette state. 



PACS numbers: 75.40.Mg, 75.10.Jm, 02.70.Lq 
I. INTRODUCTION 



The Density Matrix Technique (DMRG) has proven to 
be a very efficient method to determine the groundstate 
properties of low dimensional systemaJ. For a quantum 
chain it produces extremely accurate values for the en- 
ergy and the correlation functions. In two dimensional 
systems the calculational effort increases rapidly with the 
size of the system. The most favorable geometry is that 
of a long small strip. In practice the width of the strip 
is limited to around 8 to 10 lattice sites. Greens Func- 
tion Monte Carlo (GFMC) is not directly limited by the 
size of the system but by the efficiency of the importance 
sampling. When the system has a minus sign problem 
the statistics is ruined in the long run and accurate esti- 
mates are impossible. Many proposalsB have been made 
to alleviate or avoid the minus sign problem with varying 
success, but all of them introduce uncontrollable errors in 
the sampling. In the DMRG calculation of the wavefunc- 
tion the minus sign problem is not manifestly present. In 
all proposed cures of the minus sign problem the errors 
decrease when the guiding wavefunction approaches the 
groundstate. 

The idea of this paper is that DMRG wavefunctions are 



much better, also for larger systems, than the educated 
guesses which usually feature as guiding wave functions. 
Moreover DMRG is a general technique to construct a 
wavefunction without knowing too much about the na- 
ture of the groundstate, with the possibility to systemat- 
ically increase the accuracy. Thus DMRG wavefunctions 
would do very well when they could be used as guid- 
ing functions in the importance sampling of the GFMC. 
There is a complicating factor which prevents a straight- 
forward implementation of this idea due to the fact that 
interesting systems are so large that it is impossible to 
use a wavefunction via a look-up table. The value of the 
wavefunction in a configuration has to be calculated by 
an in-line algorithm. This has limited the guiding wave- 
functions to simple expressions which are fast to evalu- 
ate. Consequently such guiding wavefunctions are not 
an accurate representation of the true groundstate wave- 
function, in particular if the physics of the groundstate is 
not well understood. In this paper we describe a method 
to read out the DMRG wavefunction in an efficient way 
by using a special representation of the DMRG wavefunc- 
tion. 

A second problem is that a good guiding wavefunction 
alleviates the minus sign problem, but cannot remove it 
as long as it is not exact. We resolve this dilemma by ap- 
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plying the method of Stochastic Reconfiguration which 
has recently been proposed by Sorellau. The viability of 
our method is tested for the frustrated Heisenberg model. 

The behavior of the Heisenberg antiferromagnet has 
been intruiging for a long time and still is in the cen- 
ter of research. The groundstate of the antiferromag- 
netic 1-dimensional chain with nearest neighbor coupling 
is exactly known. In higher dimensions only approximate 
theories or simulation results are available. The source 
of the complexity of the groundstate are the large quan- 
tum fluctuations which counteract the tendency of classi- 
cal ordering. The unfrustrated 2-dimensional Heisenberg 
antiferromagnet orders in a Neel state and by numeri- 
cal methods the properties of this state can be analyzed 
accuratelyQ. The situation is worse when the interactions 
are competing as in a 2-dimensional square lattice with 
antiferromagnetic nearest neighbor Ji and next nearest 
neighbor J2 coupling. This spin system with continuous 
symmetry can order in 2 dimensions at zero temperature, 
but it is clear that the magnetic order is frustrated by the 
opposing tendencies of the two types of interaction. The 
ratio J2/J1 is a convenient parameter for the frustration. 
For small values the system orders antiferromagnetically 
in a Neel type arrangement, which accomodates the near- 
est neighbor interaction. For large ratios a magnetic or- 
der in alternating columns of aligned spins (columnar 
phase) will prevail; in this regime the roles of the two 
couplings are reversed: the nearest neighbor interaction 
frustrates the order imposed by the next nearest neigh- 
bor interaction. In between, for ratios of the order of 0.5, 
the frustration is maximal and it is not clear which sort 
of groundstate results. This problem has been attacked 
by various methodaJDut not yet by DMRG and only very 
recently by GFMCO. This paper addresses the issue by 
studying the spin correlations. 

A simple road to the answer is not possible since the 
behavior of the system with frustration presents some 
fundamental problems. The most severe obstacle is that 
frustration implies a sign problem which prevents the 
straightforward use of the GFMC simulation technique. 
Moreover the frustration substantially complicates the 
structure of the groundstate wavefunction. Generally 
frustration encourages the formation of local structures 
such as dimers and plaquettes which are at odds, but not 
incompatible, with long range magnetic order. These 
correlation patterns are the most interesting part of the 
intermediate phase and the main goal of this investiga- 
tion. 

Many attempts have been made to clarify the situation. 
Often simple approximations such as mean-field or spin- 
wave theory give useful information about the qualitative 
behavior of the phase diagram. A fairly sophistocated 
mean-field theory using the Schwinger boson represen- 
tation does not give an intermediate phaseQ. Given the 
complexity of the phase diagram and the subtlety of the 
effects it is not clear whether such approximate meth- 
ods can give in this case a reliable clue to the qualitative 
behavior of the system. 



Exact calculations have been performed on small sys- 
tems up to size 6 X 6 by Schulz et al.El. Although this infor- 
mation is very accurate and unbiased to possible phases, 
the extrapolation to larger systems is a long way, the 
more so in view of indications that the anticipated finite 
size behavior only applies for larger systems. Another 
drawback of these small systems is that the groundstate is 
assumed to have the full symmetry of the lattice. There- 
fore the symmetry breaking, associated with the forma- 
tion of dimers, ladders or plaquettes, which is typical for 
the intermediate state, can not be observed directly. 

More convincing are the systenaatic series expansiouj^ 
reported recently by Kotov et al.El,EI and by Singh et al.li^, 
which bear on an infinite system. They start with an 
independent dimers (plaquettes) and study the series ex- 
pansion in the coupling between the dimers (plaquettes). 
By the choice of the state, around which the perturbation 
expansion is made, the type of spatial symmetry breaking 
is fixed. These studies favor in the intermediate regime 
the dimer state over the plaquette state. Their dimer 
state has dimers organized in ladders in which the chains 
and the rungs have nearly equal strength. So the system 
breaks the translational invariance only in one direction. 
The energy differences are however small and the series 
is finite, so further investigation is useful. Our simula- 
tions yield correlations in good agreement with theirs, 
but do not confirm the picture of translational invariant 
ladders. Instead we find an additional weaker symmetry 
breaking along the ladders, such that we come closer to 
the plaquette picture. |-| 

Very recently Capriotti and SorellaS have carried out 
a GFMC simulation for J2 = 0.5Ji and have studied 
the susceptibilities for the orientational and translational 
symmetry breaking. They conclude that the groundstate 
is a plaquette state with full symmetry between the hor- 
izontal and vertical direction. 

From the purely theoretical side-the problem has been 
discussed by Sachdev and ReadO on the basis of a 
large spin expansion. From their analysis a scenario 
emerges in which the Neel phase disappears upon in- 
creasing frustration in a continuous way. Then a gapped 
spatial-inhomogeneous phase with dimerlike correlations 
appears. For even higher frustration ratios a first order 
transition takes place to the columnar phase. Although 
this scenario is qualitative, without precise location of 
the phase transition points, it definitively excludes dimer 
formation in the magnetically ordered Neel and columnar 
phase. It is remarkable that two quite different order pa- 
rameters (the magnetic order and the dimer order) disap- 
pear simultaneously and continuously on opposite sides 
of the phase transition. In this scenario, this is taken as 
an indication of some kind of duality of the two phases. 

Given all these predictions it is of utmost interest to 
further study the nature of the intermediate state. Due 
to the smallness of the differences in energy between the 
various possibilities, the energy will not be the ideal test 
for the phase diagram. Therefore we have decided to 
focus directly on the spin correlations as a function of 
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the ratio 3^1 J \- In this paper we first investigate the 2- 
dimensional frustrated Heisenberg model by construct- 
ing the DMRG wave function of the groundstate for long 
strips up to a width of 8 sites. The groundstate energy 
and the spin stiffness which are calculated, confirm the 
overal picture described above, but the results are not ac- 
curate enough to allow for a conclusive extrapolation to 
larger systems. Then we study an open 10x10 lattice by 
means of the GFMC technique using DMRG wavefunc- 
tions as guiding wavefunction for the importance sam- 
pling. The GFMC are supplcmcntedAvith Stochastic Re- 
configuration as proposed tw SorellaS as an extension of 
the Fixed Node techniqueli3. This method avoids the 
minus sign problem by replacing the walkers regularly 
by a new set of positive sign with the same statistical 
properties. The first observation is that GFMC improves 
the energy of the DMRG in a substantial and systematic 
way as can be tested in the unfrustrated model where 
sufficient information is available from different sources. 
Secondly the spin correlations become more accurate and 
less dependent on the technique used for constructing the 
DMRG wavefunction. The DMRG technique is focussed 
on the energy of the system and less on the correlations. 
The GFMC probes mostly the local correlations of the 
system as all the moves are small and correspond to local 
changes of the configurations. With these spin correla- 
tions we investigate the phase diagram for various values 
of the frustration ratio Jij J\. 

The paper begins with the definition of the model to 
avoid ambiguities. Then a short description of our im- 
plementation of the DMRG method is given. We go into 
more detail about the way how the constructed wave- 
functions can be used as guiding wavefunctions in the 
GFMC simulation. This is a delicate problem since the 
full construction of a DMRG wavefunction takes several 
hours on a workstation. Therefore we separate off the 
construction of the wavefunction and cast it in a form 
where the configurations can be obtained from each other 
by matrix operations on a vector. So the length of the 
computation of the wavefunction in a configuration scales 
with the square of the number of states included in the 
DMRG wavefunction. But even then the actual construc- 
tion of the value of the wavefunction in a given config- 
uration is so time consuming that utmost effiency must 
be reached in obtaining the wavefunction for successive 
configurations. The remaining sections are used to out- 
line the GFMC and the Stochastic Reconfiguration and 
to discuss the results. We concentrate on the correlation 
functions since we see them as most significant for the 
structure of the phases. We give first a global evalua- 
tion of the correlation function patterns for a wide set of 
frustration ratios and then focus on a number of points 
to see the dependence on the guiding wavefunction and 
to deduce the trends. The paper closes with a discussion 
and a comparison with other results in the literature. 



II. THE HAMILTONIAN 

The hamiltonian of the system refers to spins on a 
square lattice. 

W = Ji^S,-S,+J2^S, -S,. (1) 

The Si are spin \ operators and the sum is over pairs 
of nearest neigbors and over pairs of next nearest 

neighbors on a quadratic lattice. Both coupling con- 
stants Ji and Ji are supposed to be positive. J\ tries to 
align the nearest neigbor spin in an antiferromagnetic 
way and Ji tries to do the same with the next nearest 
neighbors. So the spin system is frustrated, implying an 
intrinsic minus sign in the simulations that cannot be 
gauged away by a rotation of the spin operators. 
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FIG. 1. The interaction constants Ji and J 2 



In order to prepare for the representation of the hamil- 
tonian we express the spin components in spin raising and 
lowering operators 

S, . S, = \{StS- + SrS^) + SIS]. (2) 

We will use the z component representation of the spins 
and a complete state of the spins will be represented as 

l-R) = |si, S2, • ■ ■ , sat), (3) 

where the are eigenvalues of the S] operator. The 
diagonal matrix elements of the hamiltonian are in the 
representation (^) given by 

= Jl^S.S, + J2^S,Sj. (4) 

The off-diagonal elements are between two nearby config- 
urations i?' and R. R' is the same as R except at a pair 
of nearest neighbors sites or next nearest neighbor 
sites [i, j], for which the spins Si and Sj are opposite. In 
R' the pair is turned over by the hamiltonian. Then 
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{R'\n\R) = ^Ji or {R'\n\R) = ^Ja, (5) 

depending on whether a nearest or a next nearest pair is 
flipped. 



In practice we do not solve the eigenvalues of the den- 
sity matrix in the configuration representation, but in a 
projection on a smaller basis. WhiteB has shown that 
the best way to represent the state | $) is to select the m 
eigenstates \a) with the largest eigenvalue 



III. THE DMRG PROCEDURE 

The DMRG procedure approximates the groundstate 
wavefunction by searching through |j:arious representa- 
tions in bases of a given dimension rna. Here we take the 
standard method (with two connecting sites) for granted 
and make the preparations for the extraction of the wave- 
function. 
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{Ri\p\R'i){R'M^K{Ri\a). 



(7) 




FIG. 2. The DMRG procedure with one connecting site 

The system is mapped onto a 1-dimensional chain (see 
Fig. ^ and separated into two parts: a left and right 
hand part. They are connected by one site. Each part 
is represented in a basis of at most m states. With a 
representation of all the operators in the hamiltonian in 
these bases one can find the groundstate of the system. 
We thus have several representations of the groundstate 
depending on the way in which the system is divided up 
into subsystems. The point is to see how these repre- 
sentations are connected and how they possibly can be 
improved. We take a representation for the right hand 
parts and improve those on the left. So we assume that 
for a given division we have the groundstate of the whole 
system and we want to enlarge the left hand side at the 
expense of the right hand side. The first step is to include 
the connecting site in the left hand part. This enlarges 
the basis for the left hand side from m to 2m and a se- 
lection has to be made of m basis states. This goes with 
the help of the density matrix for the left hand side as 
induced by the wavefunction for the whole system. For 
later use we write out the basic equations for the density 
matrix in the configuration representation. Let, at a cer- 
tain stage in the computation, |$) be the approximation 
to the groundstate. The configurations of the right hand 
part and the left hand part are denoted by Rr and i?;. 
Then the density matrix for the left hand part reads 



The next step is to break up the right hand part into a 
connecting site and a remainder. With the basis for this 
remainder and the newly acquired basis for the left hand 
part we can again compute the groundstate of the whole 
system as indicated in the lower part of the figure. Now 
we are in the same position as we started, with the dif- 
ference that the connecting site has moved one position 
to the right. Thus we may repeat the cycle till the right 
hand part is so small that it can exactly be represented 
by m states or less. Then we have constructed for the left 
hand part a new set of bases, all containing m states, for 
system parts of variable length. Next we reverse the roles 
of left and right and move back in order to improve the 
bases for the right hand parts with the just constructed 
bases for the left hand part. 

The process may be iterated till it converges towards 
a steady state. The great virtue of the method is that 
it is variational. In each step the energy will lower till 
it saturates. In l-dimcnsiwial system the method has 
proven to be very accurateu. So one wonders what the 
main trouble is in higher dimensions. 





{Rl\p\R[) ^Y.^Rl,Rr\^){m[.Rr 
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(a) (b) 
FIG. 3. Two 1-dimensionaI paths through the system: 
"straight" (a) and "meandering" (b). 

In Fig. ^ we have drawn 2 possible ways to map the 
system on a 1-dimensional chain. One sees that if we 
divide again the chain into a left hand part and a right 
hand part and a connecting site, quite a few sites of the 
left hand part are nearest or next nearest neighbors of 
sites of the right hand part. So the coupling between the 
two parts of the chain is not only through the connecting 
site but also through sites which are relatively far away 
from each other in the 1-dimensional path. The operators 
for the spins on these sites are not as well represented as 
those of the connecting site, which is fully represented 
by the two possible spin states. Yet the correlations be- 
tween the interacting sites count as much for the energy 
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of the system as those interacting with the connecting 
site. One may say that the further away two interacting 
sites are in the 1-dimensional chain the poorer their in- 
fluence is accounted for. This consideration explains in 
part why open systems can be calculated more accurately 
than closed systems, even in 1-dimensional systems. 

It is an open question which map of the 2-dimensional 
onto a 1-dimensional chain gives the best representation 
of the groundstate of the system. Also other divisions 
of the system than those suggested by a map on a 1- 
dimensional chain are possible and we have been exper- 
imenting with arrangements whichjxeflect better the 2- 
dimensional character of the latticetJ. They are promis- 
ing but the software for these_is not as sophisticated 
as the one developed by WhiteEfl for the 1-dimensional 
chain. We therefore have restricted our calculations to 
the two paths shown here. The second choice, the "me- 
andering" path, was motivated by the fact that it has 
the strongest correlated sites most nearby in the chain 
and this choice was indeed justified by a lower energy for 
a given dimension m of the representation than for the 
"straight" path. 

The DMRG calculations as well as the corresponding 
GFMC simulations are carried out for both paths. The 
meandering path has to be preferred over the straight 
path as the DMRG wavefunctions generally give a better 
energy value and the simulations suffer less from fluctua- 
tions. Nevertheless we have also investigated the straight 
path, since the path chosen leaves its imprints on the 
resulting correlation pattern and the paths break the 
symmetries in different ways. Both paths have an ori- 
entational preference. In open systems the translational 
symmetry is broken anyway, but the meandering path 
has in addition a staggering in the horizontal direction. 
This together with the horizontal nearest neighbor sites 
appearing in the meandering path gives a preference for 
horizontal dimerlike correlations in this path. On the 
other hand the straight path prefers the dimers in the ver- 
tical direction. Comparing the results of the two choices, 
allows us to draw further conclusions on the nature of 
the intermediate state. 



IV. EXTRACTING CONFIGURATIONS FROM 
THE DMRG WAVEFUNCTION 



It is clear that the wavefunction which results from a 
DMRG-procedure is quite involved and it is not simple 
to extract its value for a given conflguration. We assume 
now that the DMRG wavefunction has been obtained by 
some procedure and we will give below an algorithm to 
obtain efficiently the value for an arbitrary conflguration 
(see alsoEj for an alternative description). 

The first step is the construction of a set of representa- 
tions for the wavefunction in terms of two parts (without 
a connecting site in between) . Let the left hand part con- 
tain I sites and the other part N — I sites. We denote the 



TO basis states of the left hand part by the index a and 
those of the right hand part by a. The eigenstates of the 
two parts are closely linked and related as follows 



{Ri\a) = 
{Rr\a) = 



l=Y,{mi,Rr){Rr\a), 

1 (8) 



/A 



Ri 



It means that for every eigenvalue there is and eigen- 
state a for the left hand part and an a for the right hand 
density matrix. The proof of (^ follows from insertion 
in the density matrix eigenvalue equation (^. 

The second step is a relation for the groundstate wave- 
function in terms of these eigenfunctions. Generally we 
have 



(i?/,i?,|$) = J2{^i\a){Rr\(3){aP\^) 

a,f3 

while due to (||) we find 

(a^|$) = ^ (a|i?/)(/3|i?,)(i?/,i?,|$) 

Ri .Rj. 

a) = Sa.jS 



(9) 



(10) 



Thus we can represent the groundstate as 



{Rl,Rrm = 



Xl^{Ri\a)i{Rr\a)N-i- 



(11) 



For this part of the problem we have to compute and 
store the set of m eigenvalues A„ for each division /. We 
point out again that we have on the left hand side the 
wavefunction and on the right hand side representations 
for given division I, which all lead to the same wavefunc- 
tion. The last step is to see the connection between these 
representations. 

As intermediary we consider a representation of the 
wavefunction with one site si separating the spins 
si ■ ■ ■ si-i on the left hand side from s/_|_i • • ■ sn on the 
right hand side. Using the same basis as in (0) we have 



(12) 



(Sl • • • S/_i,S/, S; + i • • -SArl*) = 



We compare this representation in two ways with (|l] 
First we combine the middle site with the left hand part. 
This leads to m states which can be expressed as linear 
combinations of the states of the enlarged segment 

Y,{s, ■ ■ ■ si.,\a)^l^,isO = ^(si • • • si\a")Tl„^^,. 

a ct" 

(13) 

In fact this relation is the very essence of the DMRG 
procedure. The wave function in the larger space is pro- 
jected on the eigenstates of the the density matrix of 
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that space. Since the process of zipping back forth has 
converged there is indeed a fixed relation (|l3|). However 
when we insert (^3|) into ( [T^ ) and compare it with (|l] 
we conclude that the matrix T must be diagonal 



^a".a' — "a" ,a' \/ ^r>' ■ 



This leads to the recursion relation 



[Si - ■ ■ Sl\ 



with 



<o'(sO='/'U'(s/)/V^a 



(14) 



(15) 



(16) 



The second combination concerns the contraction of the 
middle site with the right hand part. This leads to the 
recursion relation 



[si - ■ ■ sn 



• • • SAr|a 



with 



B 



i-i 



i-i 



(17) 



(18) 



The A and B matrices are the essential ingredients of 
the calculation of the wavefunction. From (^8|) and (|l^ ) 
follows that they are related as 



<a'(«) = \A[7AF<„,(s). 



(19) 



By the recursion relations the basis states are ex- 
pressed as products of m x to matrices. The determi- 
nation of the DMRG wavefunction and the matrices A 
(or B) is part of the determination of the DMRG wave- 
function which is indeed lengthy but fortunately no part 
of the simulation. The matrices can be stored and con- 
tain the information to calculate the wavefunction for 
any configuration. The value of the wavefunction is now 
obtained as the product of matrices acting on a vector. 
Thus the calculational effort scales with to^. Using re- 
lation (|l^) one reconfirms by direct calculation that the 
wavefunction is indeed independent of the division I. 

When the simulation is in the configuration R, all the 
{Ri\a)i and the {Rr\a) n-i are calculated and stored, 
with the purpose to calculate the wavefunctions more 
efficiently for the configurations R' which are connected 
to R by the hamiltonian and which are the candidates 
for a move. The structure of of these nearby states is 
R' = si • ■ ■ ■ ■ ■ si-^ ■ ■ ■ spf {I2 > h). So we have that 
for R' the representation 



(R'm = 



E 



Aq (si • • • s;^ • • • si^\a){si, 



+1 



' SN\a) 



(20) 



holds. Now we see the advantage of having the wave- 
function stored for all the divisions. The second factor 
in (^0|) is already tabulated; the first factor involves a 
number of matrix multiplications equal to the distance 
in the chain of the two spins li and I2 till one reaches a 
tabulated function. One can use the tables for a certain 
number of moves but after a while it starts to pay off to 
make a fresh list. 



V. GREEN FUNCTION MONTE CARLO 
SIMULATIONS 

The GFMC technique employs the operator 

g = l-en (21) 
and uses the fact that the groundstate |4'o) results from 



1*0) - g^m 



e < 1, 



ne > 1 



(22) 



where in principle |<i>) may be any function which is non- 
orthogonal to the groundstate. In view of the possible 
symmetry breaking, the overlap is a point of serious con- 
cern on which we come back in the discussion. In practice 
we will use the best |<I>) that we can construct conve- 
niently by the DMRG-procedure described above. The 
closer |<I>) is to the groundstate the smaller the number 
of factors n in the product needs to be in order to find 
the groundstate. Evaluating ( p2|) in the spin represen- 
tation gives for the projection on the trial wavefunction 
the following long product 



($1*0) ^ E(^i^*^) 



R 



M 



Y[{R45\R. 



{Rom- (23) 



Here the sum is over paths R = {Rm, ■ ■ ■ Ri,Ro) which 
will be generated by a Markov process. The Markov pro- 
cess involves a transition probability T{Ri ^ Ri-i) and 
the averaging process uses a weight m{R). Its is natural 
to connect the transition probabilities to the matrix ele- 
ments of the Greens Function Q. But here comes the sign 
problem into the game: the transition probabilities have 
to be positive (and normalized). So we put the transi- 
tion rate proportional to the absolute value of the matrix 
element of the Greens Function 



T{R ^ R') = 



\{R\Gm 



Er"\{r"\s\r': 



(24) 



This implies that we have to use a sign function s(i?, R) 



s{R, R) 
and a weight factor 



{Rm') 

\{R\G\R')\ 



TO(i?)=^|(i?'|C?|i?)| 



(25) 



(26) 
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All these factors together form the matrixelement of the 
Green Function 



{R\g\R') = T{R ^ R')s{R, R')m{R') 



(27) 



If the matrix elements of the Greens Function were all 
positive, or could be made positive by a suitable trans- 
formation, we would not have to introduce the sign func- 
tion. We leave its consequences to the next section. By 
the representation ( p7| ) we can write the contribution of 
the path as a product of transition probabilities, signs 
and local weights. The transition probabilities control 
the growth of the Markov chain. The signs and weights 
constitute the weight of a path 



M{'R)=mf{RM) 



M 



\{s{R,,R,^i)m{R,^i) 



m^iRo). 
(28) 



The initial and final weight have to be chosen such that 
the weight of the paths corresponds to the expansion 
(^). For the innerproduct (<l'|4'o) we get 



(i?) = ($|i?). 



(29) 



With this final weight we have projected the groundstate 
on the trial wave. This allows us to calculate the so- 
called mixed averages. For that purpose we define the 
local estimator O 



0{R) 



mo\R) 



which yields the mixed average 



{0)r 



{'^>\0\^o) _ ERO(i?M)M(R) 



(30) 



(31) 



For operators not commuting with the hamiltonian the 
mixed average is an approximation to the groundstate 
average. Later on we will improve on it. 

In this raw form the GFMC would hardly work because 
all paths are generated with equal weight. One can do 
better by importance sampling in which one transforms 
the problem to a Greens Function with matrix elements 



{R\gm 



mR){R\Gm 
im') 



(32) 



Generally this can be seen as a similarity transformation 
on the operators and from now on everywhere operators 
with a tilde are related to their counterpart without a 
tilde as in (^). It gives only a minor change in the for- 
mulation. The transition rates are based on the matrix 
elements of G and so are the signs and weights. Thus we 
have a set of definitions like (|4|)~(|2^) with everywhere a 
tilde on top. It leads also to a change of the initial and 
final weight 



m{R) = \mR)\^' 



rhf{R) = 1. 



(33) 



By chosing these weights the formula (^0|) for the aver- 
age still applies with a weight M (R) made up as in ( [2^ ) 
with the weights and signs with a tilde. Using the tilde 
operators the local estimator ( |30| ) reads 



0{R)^Y.(R'\d\R). 



(34) 



We will speak about the various paths in terms of inde- 
pendent walkers that sample these paths. As some walk- 
ers become more important than others in the process, 
it is wise to improve the variance by branching, which 
we will discuss later with the sign problem. Before we 
embark on the discussion of the sign problem we want to 
summarize a number of aspects of the GFMC simulation 
relevant to our work. 

• The steps in the Markov process are small, only 
the ones induced by one term of the hamiltonian 
feature in a transition to a new state. This makes 
the subsequent states quite correlated. So many 
steps have to be performed before a statistical in- 
dependent configuration is reached; on the average 
a number of the order of the number of sites. 

• In every configuration the wave function for a num- 
ber of neighboring states (the ones which are reach- 
able by the Greens Function), has to be evaluated. 
This is a time consuming operation and it makes 
the simulation quite slow, because out of the pos- 
sibilities (of order N) only one is chosen and all 
the information gathered on the others is virtually 
useless. 

• The necessity to choose a small e in the Greens 
Function seems a further slow down of the method, 
but it can be avoided by the technique of contiguous 
time steps developed by Ceperley and Triveditj. In 
this method the possibility of staying in the same 
configuration (the diagonal element of the Greens 
Function) is eliminated and replaced by a waiting 
time before a move to another state is made. (For 
further pdetails in relation to the present paper we 
refer toEJ) 



• The average 

by 



can be improved by replacing it 



(35) 



of which the error with respect to the true aver- 
age is of second order in the deviation of \^) from 
l^'o). For conserved operators, such as the energy, 
this correction is not needed since the the mixed 
average gives already the correct value. 
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VI. THE SIGN PROBLEM AND ITS REMEDIES 

In the hamihonian the z component of the spin 
operator keeps the spin configuration invariant, whereas 
the X and y components change the configuration. The 
typical change is that a pair of nearest or next near- 
est neighbors is spin reversed. Inspecting the Greens 
Function it means that all changes to another configu- 
ration involve a minus sign! Thus the Greens Function 
is as far as possible from the ideal of positive matrix ele- 
ments. The diagonal terms are positive, but they always 
are positive for sufhciently small e. Importance sampling 
can remove minus signs in the transition rates, when the 
ratio of the guiding wavefunction involves also a minus 
sign. For J2 7^ no guiding wavefunction can remove 
the minus sign problem completely. In Fig. 3 we show a 
loop of two nearest neighbor spin flips followed by a flip 
in a next nearest neighbor pair, such that the starting 
configuration is restored. The product of the ratios of 
the guiding wavefunction drops out in this loop, but the 
product of the three matrix elements has a minus sign. 
So at least one of the transitions must involve a minus 
sign. 
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FIG. 4. Illustration of the sign problem in the frustrated 
Heisenberg model. The shown sequence of spin flips always 
involves a sign that can not be gauged away by a different 
choice of guiding wavefunctions 

For unfrustrated systems these loops do not exist and 
one can remove the minus sign by a transformation of 
the spin operators 



-5? 



5? 



Sf 



(36) 



which leave the commutation operators invariant. Ap- 
plying this transformation on every other spin (the white 
fields of a checkerboard) all flips involving a pair of near- 
est neighbors then give a positive matrix element for the 
Greens Function. So when J2 — the appearant sign 
problem riS| transformed away. For sufficiently small J2, 
Marshall has shown that the wave function of the sys- 
tem has only positive components (after the "Marshall" 
sign flip (|3q)). So the minus sign problem is not due to 
the wave function but to the frustration. (For the Hub- 
bard model it is the guiding wave function which must 
have minus signs due to the Pauli principle, while the 
bare transition probabilities can be taken positive). 

Due to the minus sign the weight of a long path picks 
up a arbitrary sign. Generally the weights are also grow- 
ing along a path. Thus if various paths are traced out by 
a number of independent walkers, the average over the 
paths or the walkers becomes a sum over large terms of 



both signs, or differently phrased: the average becomes 
small with respect to the variance; the signal gets lost in 
the noise. |— , 

Ceperley and AldeiO constructed a method. Fixed 
Node Monte Carlo (FNMC), which avoids the minus sign 
problem at the expense of introducing an approxima- 
tion. Their method is designed for continuum systems 
and handling fermion wavefunctions. They argued that 
the configuration space in which the wavefunction has 
a given sign, say positive, is sufficient for exploring the 
properties of the groundstate, since the other half of the 
configuration space contains identical information. Thus 
they designed a method in which the walkers remain in 
one domain of a given sign, essentially by forbidding to 
cross the nodes of the wavefunction. The approximation 
is that one has to take the nodal structure of the guid- 
ing wavefunction for granted and one cannot improve on 
that, at least not without sophistocation (nodal release). 
The method is variational in the sense that errors in the 
nodal structure always raise the groundstate energy. 

It seems trivial to take over this idea to the lattice but 
it is not. The reason is that in continuum systems one 
can make smaller steps when a walker approaches a node 
without introducing errors. In a lattice system the con- 
figuration space is discrete; so the location of the node is 
not strictly defined. The important part is that, loosely 
speaking the nodes are between configurations and one 
cannot make smaller moves than displacing a particle 
over a lattice distance or flip a pair of spins. Van Bem- 
mel et al.EZl adapted the FNMC concept to lattice sys- 
tems preserving its variational character. This extension 
to the lattice suffers from the same shortcoming as the 
method of Ceperley and Alder: the "nodal" structure 
of the guiding wavefunction is given and cannot be ini=. 
proved by the Monte Carlo process. Recently SorellaEl 
proposed a modification which overcomes this drawback. 
It is based on two ingredients. 

Sorella noticed that the following effective hamiltonian 
yields also an upper bound to the energy: 



{R\n,s\R') 



( {R\H\R') if (R\H\R') < 
\ --i{R\'H\R') if {R\n\R') > (7 > 0) 

(i?|7Yeff|i?) = (i?|7i:|i?) + (l + 7)14f(i?) 

(37) 

Here the "jSign flip" potential is the same as that of ten 
Haaf et al.li3 and given by 



V.f{R) = Y.^RMR) 



(38) 



where the subscript "na" (not-allowed) on R' restricts 
the summation to the moves for which the matrix ele- 
ment of the hamiltonian is positive (p7|). 

If the guiding wavefunction were to coincide with the 
true wavefunction, the simulation of the effective hamil- 
tonian, which is sign free by construction, yields exact 
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averages. So one may expect that good guiding wave- 
functions lead to good upperbounds for the energy. This 
upperbound increases with 7, indicating that 7 = seems 
the best choice, which is the effective hamiltonian of ten 
Haaf et al.lla. That hamihonian however is a truncated 
version of the true hamiltonian in which all the dangerous 
moves are eliminated. The sign flip potential must cor- 
rect this truncation by suppressing the probability that 
the walker will stay in a configuration with a large po- 
tential. 

The second ingredient uses the fact that the hamil- 
tonian (^7|) explores a larger phase space and therefore 
contains more information than the truncated one. Par- 
allel to the simulation of the effective hamiltonian one 
can calculate the weights for the true hamiltonian. As 
we saw in the summary of the GFMC method, forcefully 
made positive transition rates still contain the correct 
weights when supplemented with sign functions. For the 
true weights of the transition probabilities as given in 
(^), the "sign function" must be chosen as 



s{R,R') = < 



-1/7 

1 - ejRlnlR) 
I l-e(i?|7i:cff|i?) 



if {R\n\R') > 
if {R\ii\R') < 

if R^ R' 



(39) 



With these "signs" in the weights a proper average can be 
calculated, but these averages suffer from the sign prob- 
lem, the more so the smaller 7 is as one sees from (39). 
So some intermediate value of 7 has to be chosen. For- 
tunately the results are not too sensitively dependent on 
7; the value 7 = 0.5 is a good compromise and has been 
taken in our simulations. 

In any simulation some walkers obtain a large weight 
and others a small one. To lower the variance branch- 
ing is regularly applied, which means a multiplication of 
the heavily weighted walkers in favor of the removal of 
those with small weight. It is not difficult to do this in 
an unbiased way. SorellaE proposed to use the branch- 
ing much more effectively in conjunction with the signs 
defined in ( ^ ) . The average sign is an indicator of the 
usefulness of the set of walkers. Start with a set of walk- 
ers with positive sign. When the average sign becomes 
unacceptably low, the process is stopped and a reconfig- 
uration takes place. The walkers are replaced by another 
set with positive weights only, such that a number of 
measurable quantities gives the same average. The more 
observables are included the more faithful is the replace- 
ment. The construction of the equivalent set requires the 
solution of a set of linear equations. With the new set 
of walkers one continues the simulation on the basis of 
the effective hamiltonian and one keeps track of the true 
weights with signs. The reconfiguration on the basis of 
some observables gives at the same time a measure for 
these obervables. Thus measurement and reconfiguration 
go together. As the number of observables that can be in- 



cluded is limited some biases are necessarily introduced. 
Sorella showed that the error in the energy of the guiding 
wave function is easily reduced by a factor jof 10, whereas 
reduction by the FNMC of ten Haaf et al.Ej rather gives 
only a factor 2. 
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FIG. 5. The energy as function of the frustration ratio 



VII. RESULTS FOR THE DMRG 

In this section we give a brief summary of the results 
of a puEfi-DMRG-calculation. Extensive details can be 
found irJlj. The systems are strips of widths up to W = 8 
and of various lengths L. They are periodic in the small 
direction and open in the long direction. The periodicity 
enables us to study the spin stiffness. We have chosen 
open boundaries in the long direction to avoid the errors 
in the DMRG wavefunction due to periodic boundaries. 
Since we have good contrqLpf the scaling behavior in L 
we extrapolate to L — > ooEj. In the small direction we 
are restricted to W = 2, 4, 6 and 8 as odd values are 
not compatible with the antiferromagnetic character of 
the system. For wider system sizes the number of states 
which has to be taken into account exceeds the possibli- 
ties of the present workstations. Our criterion is that the 
value of the energy does not drift anymore appreciably 
upon the inclusion of more states. This does not mean 
that the wavefunction is virtually exact, since the energy 
is a rather insensitive probe for the wavefunction. For 
instance correlation functions still improve from the in- 
clusion of more states. In Fig. |^ we present the energy 
as function of the ratio J2/J1, for strip widths 4,6 and 
8 together with the best extrapolation to infinite width 
systems. The figure strongly suggests that the infinite 
system undergoes a first order phase transition around 
a value 0.6. This can be attributed to the transition to 
a columnar order (lines of opposite magnetisation). It 
is impossible to deduce more information from such an 
energy curve as other phase transitions are likely to be 
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continuous with small differences in energy between the 
phases. 

The spin stiffness can be calculated with the DMRG- 
wavefunctioUpfiDr systems which are periodic in at least 
one directiontj. 



0.8 



0.6 ^ 




FIG. 6. The stiffness ps as function of the frustration ratio. 
Finite size extrapolatjsins put the region where ps vanishes 
between 0.38 and 0.62Q 

The result of the computation is plotted in Fig. 5. One 
observes a substantial decrease of ps in the frustrated re- 
gion indicating the appearance of a magnetically disor- 
dered phase. In contrast to the energy the data do not 
allow a meaningful extrapolation to large widths. The 
lack of clear finite size scaling behavior in the regime of 
small values of W prevents to draw firm conclusions on 
the disappearence of the stiffness in the middle regime. 

For the correlation functions following from the DMRG 
wavefunction we refer totJ. 



VIII. RESULTS FOR GFMC WITH SR 

We now come to the crux of this study: the simulations 
of the system with GFMC, using the DMRG wavefunc- 
tions to guide the importance sampling. All the simula- 
tions have been carried out for 10 x 10 lattice with open 
boundaries. Standardly we have 6000 walkers and we 
run the simulations for about 10** measurements. These 
measuring points are not fully independent and the vari- 
ance is determined by chopping up the simulations into 
50-100 groups, often carried out in parallel on different 
computers. We first give an overall assessment of the cor- 
relation function pattern and then analyze some values 
of the ratio Jij J\. 

In the first series we have used the guiding wavefunc- 
tion on the basis of the meandering path Fig. ||(b), be- 
cause it gives a better energy than the straight option 
(a) . The number of basis states is m = 75, which is 



small enough to carry out the simulations with reason- 
able speed and large enough that trends begin to mani- 
fest themselves. Measurements of a number of correlation 
functions are made in conjunction with Stochastic Recon- 
figuration as described in section 7. The details of these 
calculations are given in Table |. Note that the DMRG 
guiding wavefunction gives a better energy for the me- 
andering path than for the straight path for values of 
Jij 3\ up to 0.6. From 0.7 on this difference is virtually 
absent. This undoubtly has to do with the change to 
the columnar state which can equally well be realized by 
both paths. The value of e has been chosen as a compro- 
mise: independent measurements require a large e but 
the minus sign problem requires to apply often Stochas- 
tic Reconfiguration i.e. a small e. One sees that in the 
heavily frustrated region the e must be taken small. In 
fact the more detailed calculations for J2 = 0.3Ji and 
J 2 = 0.5Ji were carried out with e = 0.01. 



In Fig. 6 and 7 we have plotted a sequence of visualiza- 
tions of the correlations. From top to bottom (zig-zag) 
they give the correlations for the values of Jij 3\. In or- 
der to highlight the differences a distinction is made be- 
tween correlations which are above average (solid lines) 
and below average (broken lines). All nearest neighbor 
spin correlations shown are negative. In all the pictures 
one sees the influence of the boundaries on the spin cor- 
relations. Only 1/4 of the lattice has been pictured, the 
other segments follow by symmetry. The upper right cor- 
ner, which corresponds to the center of the lattice, is the 
most significant for the behavior of the bulk. The overall 
trend is that spatial variations in the correlation func- 
tions occur in growing size with J2/J1. On the side of 
low J1/J2 (Neel phase) one sees dimer patterns in the 
horizontal direction, they turn over to vertical dimers 
(around J2 = 0.7Ji) and rapidly disappear in the colum- 
nar phase. This is again support for the fact that the 
columnar phase is separated from the intermediate state 
by a first order phase transition. 



Open boundary conditions have the disadvantage of 
boundary effects, which make it more difficult to dis- 
tinguish between spontaneous and induced breaking of 
the translational symmetry. On the other hand for open 
boundaries, dimers, plaquettes or any other interrup- 
tion of the translational symmetry have a natural refer- 
ence frame. The correlations are not only influenced by 
the boundaries of the system, also the guiding DMRG- 
wavefunction leaves its imprint on the results. This is 
mainly due to the fact that we have only mixed esti- 
mators for the correlation functions, which show a mix 
of the guiding wavefunction and the true wavefunction. 
The improved estimator, used in these pictures, corrects 
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for this effect to linear order in the deviation. -'^ The lad- 
der like structure in the DMRG path is reflected in a 
ladder hke pattern in the correlations as an inspection of 
the correlations in the DMRG wavefunctions (not shown 
here) reveals. But ladders are clearly also present in the 
GFMC results shown in the pictures. 

In order to eliminate the influence of the guiding wave- 
function we scrutinize some of values of J2I J\ in more de- 
tail, by inspecting how the results depend on the size of 
the basis in the DMRG wavefunction and on the choice 
of the DMRG path. Since we are mostly interested in 
the behavior of the infinite lattice, we discuss mainly the 
behavior of the correlations in and around the central pla- 
quette. So we study a sequence of DMRG wavefunctions 
for TO = 32, 75, 100, 128 and 150(200) and carry out for 
each of them extensive GFMC simulations. First we look 
to the case J2 = 0, which is easy because we know that 
it must be Neel ordered and therefore it serves as a check 
on the calculations. Then we take J2 = 0.3Ji which is 
the most difficult case since it is likely to be close to a 
phase transition. Finally we inspect Ji = 0.5Ji where 
we are fairly sure that some dimerlike phase is realized. 



A. J2 = 

For the unfrustrated Heisenberg model we have sev- 
eral checkpoints for our calculations. We can find to a 
high degree of accuracy the groundstate energy and we 
are sure that the Neel phase is homogeneous, i.e. that 
the correlations show no spatial variation other than that 
of the antiferromagnet. We have two ways of estimating 
the energy of a 10 x 10 lattice. The first method is based, 
on finite size interpolation. From DMRG calculations^ 
we have an exact value for a 4 x 4 lattice, an accurate 
value for the 6x6 lattice and a good value for the 8x8 
lattice. JThere is also the very accurate calculation of 
Sandviko for an infinitely large lattice, yielding the value 
of eo = —0.669437(5). The leading finite size correction 
goes as 1/L. Including also a 1 /L^ term we have esimated 



the value for a 10 x 10 lattice as 0.629(1) and incorpo- 
rated this value in Table ||(a). We stress that this is an 
interpolation for which the value of Sandvik is the most 
important input. 

The second method is less well founded and uses the ex- 
perience that DMRG energy estimates can be improved 
considerably by extrapolating to zero truncation error. 
When plotted as function of this truncation error the en- 
ergy is often remarkably linear. In Table Il|(b) we give for 
a series of bases m = 32, 75, 100, 128 and 150, the values 
of the truncation error and the corresponding DMRG en- 
ergy per site together with the extrapolation on the basis 
of linear behavior.^ Note that the two estimates are com- 
patible. In Table ||(b) we have also listed the values of 
the GFMC simulations for the corresponding values of 
TO. They do agree quite well with these estimates in par- 
ticular with the one based on finite size scaling. We point 
out that one would have to go very far in the number of 
states in the DMRG calculation to obtain an accuracy 
that is easily obtained with GFMC. Thus the combina- 
tion of GFMC and DMRG does really better than the 
individual components. One might wonder why there is 
still a drift to lower energy values in the GFMC simula- 
tions (which is also present in the tables to come). The 
reason is that the DMRG wavefunction is strictly zero 
outside a certain domain of configurations, because the 
truncation of the basis involves also the elimination of 
certain combinations of conserved quantities of the con- 
stituing parts. The domain of the wavefunction grows 
with the size of the basis. 

Turning now to the correlations it seems that they are 
homogeneous in the center of the lattice for J2 = 0. How- 
ever a closer inspection reveals small differences. In Table 
III we list the asymmetries in the horizontal and verti- 
cal directions of the spin correlations in and around the 
central plaquette as function of the number of states. 
If we number the spins on the lattice as Sn,m with 
1 < n,TO < 10, the central plaquette has the coordi- 
nates (5,5), (5,6), (6,5) and (6,6). We then define the 
asymmetry parameters A^; and Aj, as 



A., = 



— i(^5,4 • 85^5 + S(i 4 ■ Sg^5 + 85^6 • S5J + Sg^g • Sgj) — 5(85,5 • S5,g + S6_5 • Sg^g) 



• 87,1 



• S7J 



•Sg 



•Sg 



(40) 



^Forward walking allows to make a pure estirn|6|te of the cor- 
relations, but requires much more calculations .El 

^The value for m = 150 is not in line with the others. This 
can be explained by the fact that the construction of this 
DMRG wavefunction was slightly different from the others in 
which the basis was built up gradually. 
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So Aa; is the average value of the correlations on the 
4 horizontal bonds which are connected to the central 
plaquette minus the average of the values on the 2 hori- 
zontal bonds in the plaquette. Similarly Aj, corresponds 
to the vertical direction. The values for the asymmetry 
in Tabic [II in the vertical direction are so small that 



they have no significance. Note that the anticipated de- 
crease in /S.X is slow in DMRG and therefore also slow 
in the mixed estimator of the GFMC. The improved es- 
timator (|3^) however is truely an improvement! So one 
sees that all the observed small deviations from the ho- 
mogeneous state will disappear with the increase of the 
number of states in the basis of the DMRG wavefunction. 
(In general the accuracy of the correlations is determined 
by that of the GFMC simulations. We get as variance a 
number of the order 0.01, implying twice that value for 
the improved estimator) The vanishing of Aj, and Ay 
also prove that finite size effects are small in the center 
of the 10 X 10 lattice. From these data we may con- 
clude that the GFMC can make up for the errors in the 
DMRG wavefunction for a relative low number of basis 
states. We have not carried out a similar series for the 
straight path since this will certainly show no dimers as 
will become clear from the following cases. 



B. Ja = 0.3Ji 

This case is the most difficult to analyze since it is ex- 
pected to be close to a continuous phase transition-from 
the Neel state to a dimerlike state. As is knownE2l the 
DMRG structure of the wavefunction is not very ade- 
quate to cope with the long-range correlation in the spins 
typical for a critical point. In Table IV we have presented 
the same data as in Table [11 but now for J2 = 0.3. There 



is no pattern in the energy as function of the truncation 
error 5. The decrease of the energy as function of the size 
of the basis m is in the DMRG wavefunctions is not sat- 
urated. The GFMC simulations lead to a notably lower 
energy and they do hardly show a leveling off as function 
of the basis of the guiding wavefunction. All these points 
are indicators that the DMRG wavefunction is rather far 
from convergence and that more accurate data would re- 
quire a much larger basis. As far as the staggering in the 
correlations is concerned the values for A^^ are significant, 
also because the simulation results generally increase the 
values. Those for Aj, are not small enough to be consid- 
ered as noise. Given the fact that most authors locate 
the phase transition at higher values J2 ~ 0.4 Ji we would 
expect both A's to vanish. So either the dimerlike state 
is realized for values as low as J2 = 0.3Ji or dimer for- 
mation already starts in the Neel state. 

To get more insight in the nature of the groundstate we 
have also carried out the same set of simulations on the 
straight path (a) in Fig. |^. This guiding wavefunction 
shows virtually no formation of dimers in any direction 
as can be observed from Table ^ In spite of the fact 



that the trends indicated in the table have not come to 
convergence one may draw a few conclusions from the 
comparison of the two sets of simulations. The overal 
impression is that the meandering guiding wavefunction 
represents a groundstate of a different symmetry as com- 
pared to the straight path guiding wavefunction. The 
meandering wavefunction prefers dimers in the horizon- 
tal direction and the straight wavefunction leads to some 
dimerization in the vertical direction. The difference also 
shows up in the energy, it is not only large on the DMRG 
level but it also persists at the GFMC level. We see sim- 
ilar trends in the next case. 



C. J2 = O.5J1 



By any estimate this value of the next nearest neighbor 
coupling leads to a dimerlike state if it exists at all. No 
accurate data are available on the energy of the 10 x 10 
system to compare to our results. In Table VI we list 



the data for a set of DMRG wavefunctions with bases 
m = 32, 75, 100, 128, 150 and 200. The DMRG values of 
the energy (with exception of the value for m = 32) can 
be extrapolated to zero truncation error with the limit- 
ing value Eq = —48.4(1), which corresponds very well 
with the level in the GFMC values for larger sizes of the 
basis. This indicates again that the GFMC simulations 
can make up for the shortcoming of the DMRG wave- 
function. One would indeed have to enlarge the basis to 
m of the order of 1000 in order to achieve the value of 
the energy of the simulations which use DMRG guiding 
wavefunctions with a basis of the order of 100. 

The staggering in the correlations expressed by the 
quantities A^, for the horizontal direction and Ay for the 
vertical direction, has values that are significant. If one 
looks to the contributions of the DMRG wavefunction 
and the GFMC simulation separately, one observes that 
the overall values do agree quite well, with the tendency 
that the GFMC simulations lowers the staggerring in the 
horizontal direction and slightly increases it in the ver- 
tical direction. So we may conclude that indeed in the 
groundstate of the J2 = 0.5 Ji system, the correlations of 
the spins are not translation invariant but show a stag- 
gering. However these results neither confirm the picture 
that the dimerstate is the lowest (as suggested by Kotov 
et al.El) nor that the plaquettestate is, the groundstate (as 
concluded by Capriotti and SorellaQ). We comment on 
these discrepancies further in the discussion. 

Again it is worthwhile to compare these results with a 
simulation on the basis of the straight path (a) in Fig. ||. 
Here it is manifest that the straight path prefers to have 
the dimers in the vertical direction. Again the impression 
is that the straight path leads to a different symmetry as 
compared to the meandering path. It is not only the 
different preference in the main direction of the dimers, 
also the secondary dimerization in the perpendicular di- 
rection, notably in the meandering case, is not present 
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in the straight case. The fairly large difference in energy 
on the DMRG level becomes quite small on the GFMC 
level. 



IX. DISCUSSION 

We have presented a method to employ the DMRG 
wavefunctions as guiding wavefunctions for a GFMC 
simulation of the groundstate. Generally the combina- 
tion is much better than the two individual methods. 
The GFMC simulations considerably improve the DMRG 
wavefunction. In the intermediate regime the properties 
of the GFMC simulations depend on the guiding wave- 
function as the results for two different DMRG guiding 
wavefunctions show. 

The method has been used to observe spin correlations 
in the frustrated Heisenberg model on a square lattice. In 
this discussion we focus on the intermediate region where 
the model is most frustrated and which is the "piece de 
resistance" of the present research. We see patterns of 
strongly correlated nearest neighbor spins, to be called 
dimers. To indicate what me mean by strong and weak 
we give the values in and around the central square of the 
10 X 10 lattice, for the case J2 — 0.5Ji. In Fig. 0(c) we 
have given the values of the central square extrapolated 
to an infinite lattice. 



-0.45 -0.12 



-0.42 -0.15 



(a) (b) (c) 

FIG. 7. The correlation pattern for thiH nearest spins for 
J2 = O.5J1; (a) according to Kotov et aLQ: a dimer pattern 
in which the strength of the correlation is indicated; (b) ac- 
cording to Capriotti and Sorellaa: a plaquette state and (c) 
according to this paper: an intermediate pattern in which the 
translational invariance is broken in both directions but with 
unequal strength. The values indicated are those based on 
the meandering path and the improved estimator. 

The values are based on the improved estimator and 
it is interesting to see the trends. The horizontal strong 
correlation of -0.42 is the result of the DMRG value -0.44 
and the GFMC value -0.43, while the weak bond -0.15 is 
the result of the DMRG value -0.09 and the GFMC value 
-0.12. Thus the GFMC weakens the order parameter 
associated with the staggering. For the vertical direction 
there is hardly a change from DMRG to GFMC. One has 
to go to the next decimal to see the difference. The strong 
bond equals -0.368 and is coming from the DMRG value 
-0.375 and the GFMC value -0.371, while the improved 
weak bond of -0.271 is the resulting value of -0.275 for 
DMRG and -0.273 for GFMC. 



Before we comment on this result we discuss the in- 
fluence of the choice of the guiding wave function. We 
note that for both points Ji — 0.3Ji and J2 — 0.5Ji the 
two choices for the DMRG wavefunction give different 
results. First of all the main staggering is for the me- 
andering path (b) of Fig. ^ in the horizontal direction, 
while the straight path (a) of Fig. || prefers the dimers 
in the vertical direction. There is not much difference 
in the values of the strong and weak correlations. Sec- 
ondly the straight path shows no appreciable staggering 
in the other direction, so one may wonder whether the 
observed effect for the meandering path is real. In our 
opinion this difference has to do with the effect that the 
DMRG wavefunction "locks in" on a certain symmetry. 
The straight path yields a groundstate which is truely 
dimerlike in the sense that it is translational invariant 
in the direction perpendicular to the dimers. The me- 
andering path locks in on a different groundstate which 
holds the middle between a dimerlike and a plaquettelike 
state. The GFMC simulations cannot overcome this dif- 
ference in symmetry, likely because the two lowest states 
with different symmetry are virtually orthogonal. On the 
DMRG level there is a large difference in energy between 
the two states, favoring the meandering path strongly, on 
the GFMC level this difference has become very small. 
With this observation in mind we compare our result with 
other findings. 

The results of the series expansionsi,i and0 are shown 
111 Fig. 0(a). Their correlations organize themselves in 
spinladders. The correlations on the rungs of the ladder 
are —0.45 ± 0.5 which compares well with our strongest 
horizontal correlation and this holds also for the weak 
horizontal correlation (-0.12 vs -0.15). The most notice- 
ble difference is the value of our weak correlation in the 
vertical direction (-0.27 vs -0.36) while the strong cor- 
relation (-0.37 vs -0.36) agrees. There is no real conflict 
between our result and theirs since the symmetry they 
find is fixed by the state around which the series expan- 
sion is made. So our claim is only that our state with 
different symnietry is the lower one. In fact in the paper 
of Singh et al.t2l, it is noted that the susceptibility to a 
staggering operator in the perpendicular direction (our 
Ay) becomes very large in the dimer state for J2 = 0.5 Ji 
which we take as an indication of the nearby lower state. 
The analytical calculations inl3 andcl however do not sup- 
port the existence of the state we find. 

Neithsjj do we find support for the plaquette state 
found inla, which we have sketched in Fig. 0(b). The 
evidence of this investigation is based on the bounded- 
ness of the susceptibility for the operator which breaks 
the orientational symmetry and the divergence of the sus- 
ceptibility for the order parameter breaking translational 
invariance (corresponding to A^,). They have not sepa- 
rately investigated the values of A^, and Ay since their 
groundstate has the symmetry of the lattice and one 
would find automatically the same answer. They con- 
clude that in absence of orientational order parameter 
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and with the presence of the translational order param- 
eter the state must be plaquettehke. We beheve that 
their result is influenced by the guiding wavefunction for 
which the one-step Lanczos approximation is taken. This 
wavefunction certainly has the symmetry of the square 
and again GFMC cannot find a groundstatc with a dif- 
ferent symmetry. 

Finally we comment on the fact that we find the dimer- 
ization already for values as low as J2 = 0.3 Ji at least for 
the meandering path. As we have mentioned earlier the 
results as function of the number of states have not suf- 
ficie ritly converged to make a firm conclusion, the more 
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GFMC. Still it could be an indication that the phase 
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Thus many questions are left over, amongst others how 
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nation of DMRG and GFMC is a good tool to investigate 
these issues since they demonstrate ad oculos the corre- 
lations in the intermediate state. 
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7 

•J2 




Straight 


Meander 


e 


-Edmrg 


-Bgfmc 


-Edmrg 


Eqfmc 


0.0 


O.d 


-61.30 


-62.33(8) 


-61.84 


-62.54(4) 


0.1 


0.06 


-57.96 




-58.53 


-59.25(2) 


0.2 


0.04 


-54.75 


-56.08(11) 


-55.48 


-56.22(4) 


0.3 


0.02 


-51.75 


-53.17(4) 


-52.50 


-53.38(3) 


0.4 


0.02 


-49.00 


-50.51(8) 


-49.92 


-50.60(5) 


0.5 


0.014 


-46.68 


-47.76(6) 


-47.78 


-48.34(4) 


0.6 


0.015 


-45.41 




-46.03 


-46.40(3) 


0.7 


0.015 


-45.67 




-45.60 


-46.00(2) 


0.8 


0.02 


-49.16 




-49.13 


-49.60(9) 


0.9 


0.02 


-53.61 




-53.70 


-54.52(2) 


1.0 


0.02 


-58.46 


-59.71(9) 


-58.64 


-59.80(8) 



TABLE I. For each degree of frustration the imaginary time interval e, the energy of the guiding state J^dmrg and that of 
the GFMC state Eqfmc are listed. 



L 


eo{L X L) 


4 


-0.5740 


6 


-0.6031 


8 


-0.6188 


10 


-0.629(1) 


oo 


-0.669437(5) 



# states 


trunc. error 


eo (DMRG) 


eo (GFMC) 


32 


21.2 xlO"-"" 


-0.6084 


-0.6192(1) 


75 


12.0 xlO"^ 


-0.6184 


-0.6254(5) 


100 


10.5 xlO"-"^ 


-0.6201 


-0.625(2) 


128 


8.7 xlO-'^ 


-0.6214 


-0.6269(6) 


150 


9.6 xlO-5 


-0.6231 


-0.6277(5) 


2N 





-0.631(3) 





(a) (b) 
TABLE II. Interpolation (a) and extrapolation (b) estimates of the energy per site of a 10 x 10 lattice 



# states 


Ax 


A. 


m 


DMRG 


GFMC 


Improved 


DMRG 


GFMC 


Improved 


32 


0.14373 


0.09981 


0.05589 


-0.00060 


0.00078 


0.00216 


75 


0.07291 


0.05668 


0.04045 


0.00081 


0.00601 


0.01121 


100 


0.06432 


0.04255 


0.03088 


0.00030 


0.00173 


0.00316 


128 


0.05619 


0.03734 


0.01849 


0.00091 


-0.00040 


-0.00173 


150 


0.05044 


0.03612 


0.02221 


0.00079 


0.00261 


0.00442 



TABLE III. Values for the asymmetry in the center for J2 = 0. As discussed in the text the error in the improved estimator 
values is of the order 0.02, which means that for m = 128 and higher the values are statistically indistinguishable from zero. 



15 



^ states 
m 


DMRG 


GFMC 


(5* 10^ 


^'dmrg 




A^ 


Egfmg 


A:, 


A, 


32 


19.0 


-51.609 


0.27784 


0.00295 


-52.81(43) 


0.363 


-0.009 


75 


10.6 


-52.581 


0.15462 


0.00616 


-53.29(05) 


0.207 


0.011 


100 


9.4 


-52.707 


0.14709 


0.00943 


-53.32(33) 


0.145 


0.009 


128 


10.6 


-52.821 


0.13042 


0.00577 


-54.01(04) 


0.254 


0.063 


150 


10.4 


-52.888 


0.12564 


0.00737 


-54.10(12) 


0.236 


0.103 



TABLE IV. Energies and asymmetries for the case J2 = 0.3Ji as function of the number of basis states m. 5 is the 
truncation error. The asymmetries A^; and Ay for the GFMC simulations are calculated with the improved estimator. The 
guiding wavefunction is obtained from the meandering path (b) in Fig. ^. The statistical error in Aa; and A^ is of the order 
0.02 



# states 

TO 


DMRG 


GFMC 


(5* 10^ 


£^DMRG 


A^ 


Ay 


-Egfmc 


A:, 


\ 


32 


30.0 


-50.672 


0.00032 


0.01657 


-52.15(11) 


0.061 


0.047 


75 


18.9 


-51.733 


-0.00295 


0.00426 


-53.21(10) 


-0.030 


0.036 


100 


19.9 


-52.066 


0.00349 


0.00492 


-53.84(72) 


0.061 


0.079 


128 


24.6 


-52.302 


0.00139 


0.00791 


-53.50(19) 


0.079 


0.027 


150 


25.7 


-52.455 


0.00222 


0.00780 


-53.52(10) 


0.022 


0.065 



TABLE V. Comparison of the energies and the values for the asymmetry in the center for the DMRG wavefunction based 
on the first (straight) path (a) in Fig. H and the associated GFMC simulation; J2 = 0.3Ji 
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# states 
m 


DMRG 


GFMC 


6*10"" 


-Bdmrg 


A:, 


A^ 


-Bgfmc 


A^ 


A, 


32 


11.8 


-47.116 


0.43245 


0.14667 


-47.55(29) 


0.295 


0.065 


75 


17.4 


-47.771 


0.38954 


0.13059 


-48.22(04) 


0.339 


0.070 


100 


12.4 


-47.924 


0.39364 


0.07877 


-48.37(22) 


0.310 


0.110 


128 


8.4 


-48.014 


0.37317 


0.08246 


-48.32(05) 


0.336 


0.139 


150 


8.3 


-48.088 


0.35819 


0.07983 


-48.33(12) 


0.324 


0.112 


200 


7.6 


-48.153 


0.34590 


0.09973 


-48.43(05) 


0.272 


0.094 



TABLE VI. Energies and asymmetries for J2 = 0.5Ji with guiding wavefunction based on the meandering path (b) in Fig. 



# states 
m 


DMRG 


GFMC 


5*10^ 


i?DMRG 


A^ 


Ay 


i?GFMC 


A^ 




32 


69.4 


-45.756 


0.00172 


0.24701 


-47.45(08) 


0.074 


0.185 


75 


26.2 


-46.718 


0.00171 


0.34950 


-47.81(25) 


-0.025 


0.302 


100 


21.2 


-46.993 


0.00063 


0.33131 


-48.16(06) 


-0.003 


0.350 


128 


24.6 


-47.231 


-0.00029 


0.32994 


-48.31(08) 


0.013 


0.291 


150 


25.7 


-47.379 


0.00215 


0.32458 


-48.33(06) 


-0.026 


0.257 



TABLE VII. Same as Table ^ but now for the "straight" path Fig (a) | 
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FIG. 8. The relative correlation strengths on 10 x 10 lattice. All other nearest neighbour correlations can be obtained by 
reflection these picture in the two dashed lines. The DMRG guiding state follows the meandering sequence of Fig. ^(b). More 
explanation is given in the text. Reading zig zag from top left to bottom right, the values for J2 are J2 = 0, . . . , 0.5 in steps of 
0.1. 
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FIG. 9. The continuation of figure pi the relative correlation strengths on 10 x 10 lattice. J2 — 0.5, . . . , 1.0 in steps of 0.1. 
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